Data description

Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.
In this file we have compared the 3 models HMSC,GJAM, JSDM on 21 simulated Dataset.\ We have compared the ability to recover the true interactions used for modelling this datasets, by comparing estimated correlation matrix and matrix representing true interactions.
There is a separate file, with the functions needed to fit each model and study the convergence and more parameters are presented.
## Environment filtering 5 species

JSDM

data<-sim_data$EnvEvenSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

Y_cor<-cor(data$Y)

to_prec<-function(m){
  n<-dim(m)[1]
  Tau_n<-matrix(nrow=n, ncol=n)
  for (j in 1:n) {
   for (k in 1:n){
    Tau_n[j, k] <-  -m[j, k]/sqrt((m[j,j]*m[k,k]))
   }
  }
  return(Tau_n)
}

#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
cols = colorRampPalette(c("dark blue","white","red"))
col2 <- colorRampPalette(c("#4393C3", "#2166AC", "#053061",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#67001F", "#B2182B", "#D6604D", "#F4A582"))


corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")

j_metric_e5<-metrics_jsdm(me5)
#cat("Success rate ")
#j_metric_e5$success_env
 cat(sprintf("Success rate: %s\n", metrics_jsdm(me5)$success_env))
## Success rate: 0.9
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =as.matrix(me5$mean$Rho*(!me5$overlap0$Rho))) 
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =as.matrix(me5$mean$Rho*(!me5$overlap0$Rho))) 

GJAM

############################################################Functions
makeSymm <- function(m) {
  m[upper.tri(m)] <- t(m)[upper.tri(m)]
  return(m)
}


expandSigma_rmd <- function(sigma, S){

  ss <- diag(S)
  ss[lower.tri(ss,diag=T)] <- sigma
  ss[upper.tri(ss)] <- t(ss)[upper.tri(ss)]
  ss
}



convert_to_m<-function(ar){
  d <-floor((sqrt(length(ar)*8+1)-1)/2)
  C <- matrix(0,d,d)
  i.lwr <- which(lower.tri(C, diag = TRUE), arr.ind=TRUE)
  C[i.lwr] <- ar
  C<-makeSymm(C)
  return(t(C))
}


metrics_gjam<-function(Rho, comp=NULL,fac=NULL,only_env=T){
    if(!is.null(comp) & !is.null(fac)){
    fac_comp <- fac-comp
    TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
    for(i in 1:(nrow(Rho)-1)){
      for(j in (i+1):ncol(Rho)){
        if(fac_comp[i,j]>0 & Rho[i,j] >0)
          TP_fac=TP_fac+1
        if(fac_comp[i,j]>0 & Rho[i,j] == 0)
          FN_fac=FN_fac+1
        if(fac_comp[i,j]==0 & Rho[i,j] == 0)
           TN=TN+1
        if(fac_comp[i,j]==0 & Rho[i,j] != 0)
           FP=FP+1
        if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
          TP_comp=TP_comp+1
        if(fac_comp[i,j]<0 & Rho[i,j] == 0)
          FN_comp=FN_comp+1
        if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
          wrong=wrong+1
      }
    }
    
  }else{
    if(!is.null(comp)){
    TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
    TP_fac <- FN_fac<-NULL
    for(i in 1:(nrow(Rho)-1)){
      for(j in (i+1):ncol(Rho)){
        if(comp[i,j]>0 & Rho[i,j] <0)
          TP_comp=TP_comp+1
        if(comp[i,j]>0 & Rho[i,j] >0)
          wrong=wrong+1
        if(comp[i,j]>0 & Rho[i,j] == 0)
          FN_comp=FN_comp+1
        if(comp[i,j]==0 & Rho[i,j] == 0)
          TN=TN+1
        if(comp[i,j]==0 & Rho[i,j] != 0)
          FP=FP+1
       
      }
    }
    }
  if(!is.null(fac)){
    
            TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
            TP_comp <- FN_comp<-NULL
            for(i in 1:(nrow(Rho)-1)){
              for(j in (i+1):ncol(Rho)){
                   if(fac[i,j]>0 & Rho[i,j] >0)
                      TP_fac=TP_fac+1
                   if(fac[i,j]>0 & Rho[i,j] <0)
                      wrong = wrong+1
                   if(fac[i,j]>0 & Rho[i,j] == 0)
                      FN_fac=FN_fac+1
                   if(fac[i,j]==0 & Rho[i,j] == 0)
                      TN=TN+1
                   if(fac[i,j]==0 & Rho[i,j] != 0)
                      FP=FP+1
              }
            }
  }
  }
  if(only_env){
    
      FP <- TN <- 0
      TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
         for(i in 1:(nrow(Rho)-1)){
            for(j in (i+1):nrow(Rho)){
                if(Rho[i,j] != 0)
                      FP=FP+1
                if(Rho[i,j] == 0)
                      TN=TN+1
              }
            }
      
  }
  success_env<-TN/(TN+FP)
  
  if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
  
  if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
  
  list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
       "wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
  }


fit_gjam<-function(data, it=2500,burn=500 , name="./gjam_models/temp.rda",interact=diag(ncol(data$Y))){
  #setup parameters
  data <- list(
    Y = subset(data, select = -env),
    X = cbind(1, scale(poly(data$env, 2))),
    covx = cov(cbind(1, scale(poly(data$env, 2)))),
    K = 3,
    J = ncol(data) - 1,
    n = nrow(data),
    I = diag(ncol(data) - 1),
    df = ncol(data)
  )
  xdata<-as.data.frame(data$X[,-1])
  colnames(xdata)<- c("env","env2")
  ydata<-as.data.frame(data$Y)
  ###formula
  formula<-as.formula( ~env+ env2)
  ml   <- list(ng = it, burnin = burn, typeNames = 'PA')
  ####fit
  
  mod_gjam1  <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
  save(mod_gjam1, file = name)
  #summary(mod_gjam1)
  hist(effectiveSize(mod_gjam1$chains$sgibbs), main="ess(beta)",lwd=2,col=gray(.6))
  hist(effectiveSize(mod_gjam1$chains$bgibbs), main="ess(beta)",lwd=2,col=gray(.6))
  
  # postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
  # postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
  postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
  postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
  
  pH<-convert_to_m(postH)
  pL<-convert_to_m(postL)
  
  R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
  
  
  sgibbs<-mod_gjam1$chains$sgibbs
  tau<-array(NA,dim=c(data$J,data$J,1000))
  for(j in 1:1000){
    ss <- expandSigma_rmd(sgibbs[j,], S = data$J)
    si <- solve(ss)
    tau[,,j] <- -cov2cor(si)
  }
  #  sgibbs<-mod_gjam1$chains$sgibbs
  # tau<-array(NA,dim=c(data$J,data$J,it-burn))
  # for(j in (burn:it)){
  #   ss <- expandSigma_rmd(sgibbs[j,], S = 5)
  #   si <- solve(ss)
  #   tau[,,j-burn] <- -cov2cor(si)
  # }

  tau_mean<-apply(tau,c(1,2), mean)
  tau_HI<-apply(tau,c(1,2),quantile,0.95)
  tau_LO<-apply(tau,c(1,2),quantile,0.05)
  
  Tau_sign<-tau_mean*(!(tau_HI>0 & tau_LO<0))
  
  # par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
  # corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200),   type = "lower")
  # title("Correlation cor(Y)")
  # corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  # title("R")
  # corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  # title("E matrix")
  # corrplot(Tau_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  # title("Tau signif")
  # corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  # title("R signif")
  # corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  # title("True interactions")
   R_sign
}


###this function only loads the model and return the R and T significant
load_gjam<-function(data,it=2500,burn=500,name="./gjam_models/temp.rda",interact=diag(ncol(data$Y)),autoC=F){
  #setup parameters
  data <- list(
    Y = subset(data, select = -env),
    X = cbind(1, scale(poly(data$env, 2))),
    covx = cov(cbind(1, scale(poly(data$env, 2)))),
    K = 3,
    J = ncol(data) - 1,
    n = nrow(data),
    I = diag(ncol(data) - 1),
    df = ncol(data)
  )

  
  gj_mod<-load_object(name)
  #summary(gj_mod)
  S<-ncol(data$Y)
  ind <-which(array(diag(S)[lower.tri(diag(S),diag=T)])==1)
  par(mfrow=c(1,3),oma = c(1, 1, 1, 1))
  hist(effectiveSize(gj_mod$chains$sgibbs[,-ind]), main="sigma non diagonal elements",lwd=2,col=gray(.6))
  hist(effectiveSize(gj_mod$chains$sgibbs[,ind]), main="sigma diagonal elements",lwd=2,col=gray(.6))
  hist(effectiveSize(gj_mod$chains$bgibbs), main="ess(beta)",lwd=2,col=gray(.6))
  
   if(autoC){
      plot(acfplot(mcmc(gj_mod$chains$bgibbs)))
      plot(acfplot(mcmc(gj_mod$chains$sgibbs)))
  }

  
  postH<-apply(gj_mod$chains$sgibbs[burn:dim(gj_mod$chains$sgibbs)[1],], 2, quantile,0.95)
  postL<-apply(gj_mod$chains$sgibbs[burn:dim(gj_mod$chains$sgibbs)[1],], 2, quantile,0.05)
  pH<-convert_to_m(postH)
  pL<-convert_to_m(postL)
  R_sign<-cov2cor(gj_mod$parameters$sigMu)*(!(pH>0 & pL<0))
  
  
  sgibbs<-gj_mod$chains$sgibbs
  tau<-array(NA,dim=c(data$J,data$J,it-burn))
  for(j in (burn:it)){
    ss <- expandSigma_rmd(sgibbs[j,], S = data$J)
    si <- solve(ss)
    tau[,,j-burn] <- -cov2cor(si)
  }
  
  tau_mean<-apply(tau,c(1,2), mean)
  tau_HI<-apply(tau,c(1,2),quantile,0.95)
  tau_LO<-apply(tau,c(1,2),quantile,0.05)
  
  Tau_sign<-tau_mean*(!(tau_HI>0 & tau_LO<0))
  
  

  par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
  corrplot(cor(data$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Correlation cor(Y)")
  corrplot(gj_mod$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R")
  corrplot(gj_mod$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("E matrix")
  corrplot(Tau_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Tau signif")
  corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R signif")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
  
  return(list(Rho_sign=R_sign,Tau_sign=Tau_sign))
}

######################################################################
##to setup chains parameters
data<-sim_data$EnvEvenSp5
#Rho_sign<-fit_gjam(data,2500,500,name="./gjam_models/gjam5env.rda")
mod_gjam<-load_gjam(data,name="./gjam_models/gjam5env.rda", autoC=T)

g_metric_e5<-metrics_gjam(mod_gjam$Rho_sign)
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_gjam(mod_gjam$Rho_sign)$success_env))
## Success rate: 0.5
mod_list_Rho<-list.append(mod_list_Rho,gjam=mod_gjam$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=mod_gjam$Tau)

HMSC

#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="Fit",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
  if (label=="Fit"){
    Y_data = subset(data, select = -env)
    ns<- ncol(Y_data)
    np <- nrow(Y_data)
    X<-scale(poly(data$env[1:np], 2))
    colnames(X)<-c("env","env2")
    studyDesign = data.frame(sample = as.factor(1:np))
    rL = HmscRandomLevel(units = studyDesign$sample)
    m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
             studyDesign = studyDesign, ranLevels = list(sample = rL))
    m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
    save(m, file = name)
    return(m)
  }
  if (label=="Load"){
    return(load_object(name))
  }
}


data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")

Convergence:

hm_conv<-function(mod, autoC=F){
  codaList = convertToCodaObject(mod)

  #convergence histograms
  hist(effectiveSize(codaList$Beta), main="ess(beta)",lwd=2,col=gray(.6))
  hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(beta)")
  if(autoC){
      plot(acfplot(codaList$Beta))
      plot(acfplot(codaList$Omega[[1]]))
  }
  hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)",lwd=2,col=gray(.6))
  hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf,lwd=2,col=gray(.6),      main="psrf(omega)")
}

hm_conv(mod=hm_mod,autoC=TRUE )

Study of interactions

metrics_hmsc<-function(Rho, comp=NULL,fac=NULL,only_env=T){
  #{
    if(!is.null(comp) & !is.null(fac)){
    fac_comp <- fac-comp
    TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
    for(i in 1:(nrow(Rho)-1)){
      for(j in (i+1):ncol(Rho)){
        if(fac_comp[i,j]>0 & Rho[i,j] >0)
          TP_fac=TP_fac+1
        if(fac_comp[i,j]>0 & Rho[i,j] == 0)
          FN_fac=FN_fac+1
        if(fac_comp[i,j]==0 & Rho[i,j] == 0)
           TN=TN+1
        if(fac_comp[i,j]==0 & Rho[i,j] != 0)
           FP=FP+1
        if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
          TP_comp=TP_comp+1
        if(fac_comp[i,j]<0 & Rho[i,j] == 0)
          FN_comp=FN_comp+1
        if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
          wrong=wrong+1
      }
    }
    
  }else{
    if(!is.null(comp)){
    TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
    TP_fac <- FN_fac<-NULL
    for(i in 1:(nrow(Rho)-1)){
      for(j in (i+1):ncol(Rho)){
        if(comp[i,j]>0 & Rho[i,j] <0)
          TP_comp=TP_comp+1
        if(comp[i,j]>0 & Rho[i,j] >0)
          wrong=wrong+1
        if(comp[i,j]>0 & Rho[i,j] == 0)
          FN_comp=FN_comp+1
        if(comp[i,j]==0 & Rho[i,j] == 0)
          TN=TN+1
        if(comp[i,j]==0 & Rho[i,j] != 0)
          FP=FP+1
      }
    }
  #  }
  #  }
  }
  if(!is.null(fac)){
    
            TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
            TP_comp <- FN_comp<-NULL
            for(i in 1:(nrow(Rho)-1)){
              for(j in (i+1):ncol(Rho)){
                   if(fac[i,j]>0 & Rho[i,j] >0)
                      TP_fac=TP_fac+1
                   if(fac[i,j]>0 & Rho[i,j] <0)
                      wrong = wrong+1
                   if(fac[i,j]>0 & Rho[i,j] == 0)
                      FN_fac=FN_fac+1
                   if(fac[i,j]==0 & Rho[i,j] == 0)
                      TN=TN+1
                   if(fac[i,j]==0 & Rho[i,j] != 0)
                      FP=FP+1
              }
            }
  }
  }  
  
  if(only_env){
    
      FP <- TN <- 0
      TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
         for(i in 1:(nrow(Rho)-1)){
            for(j in (i+1):nrow(Rho)){
                if(Rho[i,j] != 0)
                      FP=FP+1
                if(Rho[i,j] == 0)
                      TN=TN+1
              }
            }
      
  }
  success_env<-TN/(TN+FP)
  
  if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
  
  if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
  
  list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
       "wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
}



hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
  getOmega = function(a,r=1)
  return(crossprod(a$Lambda[[r]]))
  ns<-mod$ns
  postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
  postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
  
  postOmega<-abind(postOmega1,postOmega2,along=3)
  postOmegaMean = apply(postOmega,c(1,2),mean)
  postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
  postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
  
  postR<-array(dim=c(ns,ns,nchains*nsamples))
  postT<-array(dim=c(ns,ns,nchains*nsamples))
  for(i in 1:dim(postOmega)[3]){
  postR[,,i]<-stats::cov2cor(postOmega[,,i])
  postT[,,i]<- -stats::cov2cor(solve(postOmega[,,i]+diag(ns)))
  }
  
  postRMean = apply(postR,c(1,2),mean)
  postRUp=apply(postR,c(1,2),quantile,0.95)
  postRLo=apply(postR,c(1,2),quantile,0.05)
  
  postTMean = apply(postT,c(1,2),mean)
  postTUp=apply(postT,c(1,2),quantile,0.95)
  postTLo=apply(postT,c(1,2),quantile,0.05)
  
  Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
  Toplot_T<-postTMean*(!(postTUp>0 & postTLo<0))
  
  # Omegacor<- computeAssociations(m)
  # supportLevel<- 0.95
  # toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
  # corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
  
  par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
  corrplot(cor(hm_mod$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Correlation cor(Y)")
  corrplot(postRMean, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R")
  corrplot(Toplot_R, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Plot only non zero value")
  corrplot(Toplot_T, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Partial correlation matrix")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
  
  return(list(Rho_sign=Toplot_R,Tau_sign=Toplot_T)) 
}


hm_mod_R<-hm_inter(hm_mod,nchains=2,nsamples = 3000)

h_metric_e5<-metrics_hmsc(hm_mod_R$Rho_sign)
# cat("Success rate ")
# h_metric_e5$success_env
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_hmsc(hm_mod_R$Rho_sign)$success_env))
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
 par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
  corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R  JSDM ")
  corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R GJAM")
  corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R HMSC")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(5))

ALL4<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
 par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
  #corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  #title("Partial correlation  JSDM ")
  corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Partial correlation GJAM")
  corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Partial correlation HMSC")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(5))

Environmental filtering 10 species

JSDM

me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)
## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
jsdm_conv(me10)
## Maximum Rhat value for Beta: 1.5602
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.7953

me10$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
j_metric_e10<-metrics_jsdm(me10)
cat("Success rate ")
## Success rate
j_metric_e10$success_env
## [1] 0.9333333
data<-sim_data$EnvEvenSp10
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
  par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
  corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Correlation cor(Y)")
  corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("EnvRho")
  corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("EnvRho signif")
  corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Rho")
  corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Rho signif")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}

#plot_cor_jsdm(me10,data$Y)


mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =me10$mean$Rho*(!me10$overlap0$Rho)) 
 
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =me10$mean$Rho*(!me10$overlap0$Rho)) 

Gjam

data<-sim_data$EnvEvenSp10

#Rho_sign<-fit_gjam(data,3000,500,"./gjam_models/gjam10env.rda")
gjam_mod<-load_gjam(data,3000,500,name="./gjam_models/gjam10env.rda")

#gje10<-load_object("./gjam_models/gjam10env.rda")

g_metric_e10<-metrics_gjam(gjam_mod$Rho_sign)
cat("Success rate ")
## Success rate
g_metric_e10$success_env
## [1] 0.3111111
#to check posterior density of s in Sigma 
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$EnvEvenSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10env.rda" )

hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod)

h_metric_e10<-metrics_hmsc(hm_mod_R$Rho_sign)
# cat("Success rate ")
# h_metric_e10$success_env
# h_metric_e10$success_fac

cat(" ")
cat(sprintf("Success rate: %s\n", h_metric_e10$success_env))
## Success rate: 1
cat(sprintf("Success rate: %s\n", h_metric_e10$success_fac))


mod_list_Rho<-list.append(mod_list_Rho, hmsc=hm_mod_R$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, hmsc=hm_mod_R$Tau_sign) 

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(10))

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(10))

Environmental filtering 20 species

JSDM

me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)
## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
me20$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
jsdm_conv(me20)
## Maximum Rhat value for Beta: 1.2174
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3236

j_metric_e20<-metrics_jsdm(me20)
cat("Success rate ")
## Success rate
j_metric_e20$success_env
## [1] 0.9894737
data<-sim_data$EnvEvenSp20
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################

plot_cor_jsdm(me20,data$Y)

mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =me20$mean$Rho*(!me20$overlap0$Rho)) 
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =me20$mean$Rho*(!me20$overlap0$Rho)) 

Gjam

data<-sim_data$EnvEvenSp20

#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")

gjam_mod<-load_gjam(data,5000,500,name="./gjam_models/gjam20env.rda")

#gje20<-load_object("./gjam_models/gjam20env.rda")

g_metric_e20<-metrics_gjam(gjam_mod$Rho_sign)
cat("Success rate ")
## Success rate
#g_metric_e20$success_env

cat(" ")
cat(sprintf("Success rate env: %s\n", round(metrics_gjam(gjam_mod$Rho_sign)$success_env,3)))
## Success rate env: 0.358
#to check posterior density of s in Sigma 
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

Gjam dimension reduction

# 
# data<-sim_data$EnvEvenSp20
# 
# #fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
# 
# #load_gjam(data,name="./gjam_models/gjam20env.rda")
# #gje20<-load_object("./gjam_models/gjam20env.rda")
# 
# 
# #to check posterior density of s in Sigma 
# #gje20<-load_object("./gjam_models/gjam20env.rda")
# #plot(density(gje20$chains$sgibbs[,4]))
#  data <- list(
#     Y = subset(data, select = -env),
#     X = cbind(1, scale(poly(data$env, 2))),
#     covx = cov(cbind(1, scale(poly(data$env, 2)))),
#     K = 3,
#     J = ncol(data) - 1,
#     n = nrow(data),
#     I = diag(ncol(data) - 1),
#     df = ncol(data)
#   )
# xdata<-as.data.frame(data$X[,-1])
# colnames(xdata)<- c("env","env2")
# ydata<-as.data.frame(data$Y)
# ###formula
# rl   <- list(r = 8, N = 20)
# formula<-as.formula( ~env+ env2)
# ml  <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
# ####fit
# 
# 
# mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
# save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
# summary(mod_gjam1)
# 
# Tau <- solve(mod_gjam1$parameters$sigMu)
# Tau_n = to_prec(Tau)
# 
# #postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
# #postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
# 
# 
# #pH<-convert_to_m(postH)
# #pL<-convert_to_m(postL)
# 
# #R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
# 
# g_metric_e20_dim_red<-metrics_gjam(Rho_sign)
# cat("Success rate ")
# g_metric_e20_dim_red$success_env
# 
# par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
# corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("Correlation cor(Y)")
# corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("R")
# corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("E matrix")
# # corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# # title("Tau")
# # corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# # title("R signif")
# corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("True interactions")
# 
data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod)

h_metric_e20<-metrics_hmsc(hm_mod_R$Rho_sign)
cat("Success score ")
## Success score
h_metric_e20$success_env
## [1] 1
cat(" ")
cat(sprintf("Success rate non-interacting: %s\n", round(metrics_hmsc(hm_mod_R$Rho_sign)$success_env,3)))
## Success rate non-interacting: 1
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(20))

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(20))

Environmental filtering + Facilitation dense 5 species

JSDM

mfd5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mfd5)
## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
#mfd5$Rhat
jsdm_conv(mfd5)
## Maximum Rhat value for Beta: 1.3081
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1421

mfd5$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
j_metric_facDense5<-metrics_jsdm(mfd5,fac = fac_inter[[4]],only_env = F)
cat("Success rate ")
## Success rate
j_metric_facDense5$success_env
## [1] 1
j_metric_facDense5$success_fac
## [1] 0
data<-sim_data$FacDenseSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

#plot_cor_jsdm(mfd5,data$Y,fac_inter[[4]])

mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =mfd5$mean$Rho*(!mfd5$overlap0$Rho)) 
 mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =mfd5$mean$Rho*(!mfd5$overlap0$Rho)) 

Gjam

data<-sim_data$FacDenseSp5
#Rho_sign<-fit_gjam(data,10000,2000,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
gjam_mod<-load_gjam(data,10000,2000,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])

#gjfd5<-load_object("./gjam_models/gjam5f.rda")

g_metric_facDense5<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[4]], only_env = F)
#cat("Success rate ")
#g_metric_facDense5$success_env
#g_metric_facDense5$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0.5
mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign)

HMSC

data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[4]])
## $Rho_sign
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    0    0    0    0
## [2,]    0    1    0    0    0
## [3,]    0    0    1    0    0
## [4,]    0    0    0    1    0
## [5,]    0    0    0    0    1
## 
## $Tau_sign
##      [,1] [,2] [,3] [,4] [,5]
## [1,]   -1    0    0    0    0
## [2,]    0   -1    0    0    0
## [3,]    0    0   -1    0    0
## [4,]    0    0    0   -1    0
## [5,]    0    0    0    0   -1
hm_mod_R<-hm_inter(hm_mod)

h_metric_facDense5<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[4]],only_env = F)
#cat("Success rate ")
# h_metric_facDense5$success_env
# h_metric_facDense5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense5$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[4]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[4]])

Environmental filtering + Facilitation dense 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2866
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4808

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 0.925
## [1] 0
data<-sim_data$FacDenseSp10

#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
gjam_mod<- load_gjam(data,5000,500,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])

#gjfd5<-load_object("./gjam_models/gjam10fd.rda")

g_metric_facDense10<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[5]], only_env = F)
# cat("Success rate ")
# g_metric_facDense10$success_env
# g_metric_facDense10$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense10$success_env,3)))
## Success rate for non-interacting: 0.625
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense10$success_fac,3)))
## Success rate for facilitation: 0.8
mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign)

HMSC

## Success rate for non-interacting: 1
## Success rate for facilitation: 0

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[5]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[5]])

Environmental filtering + Facilitation dense 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.6132
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.508

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9833333
## [1] 0

Gjam

data<-sim_data$FacDenseSp20

#Rho_sign<-fit_gjam(data,10000,1500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])
gjam_mod<-load_gjam(data,10000,1500,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])

#gjfd5<-load_object("./gjam_models/gjam20fd.rda")

g_metric_facDense20<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[6]], only_env = F)
# cat("Success rate ")
# g_metric_facDense20$success_env
# g_metric_facDense20$success_fac
# 
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 0.406
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0.5
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[6]])

h_metric_facDense20<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[6]],only_env = F)
# cat("Success rate ")
# h_metric_facDense20$success_env
# h_metric_facDense20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense20$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[6]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[6]])

Environmental filtering + Facilitation sparse 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.373
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4048

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.889
## Success rate for facilitation: 1

Gjam

data<-sim_data$FacSparseSp5

#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam5fs.rda",interact=fac_inter[[7]])
gjam_mod<-load_gjam(data,2500,500,name="./gjam_models/gjam5fs.rda", interact=fac_inter[[7]])

#gjfd5<-load_object("./gjam_models/gjam5fs.rda")

g_metric_facSparse5<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[7]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse5$success_env
# g_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacSparseSp5

hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fs.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[7]])
h_metric_facSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse5$success_env
# h_metric_facSparse5$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse5$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[7]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[7]])

Environmental filtering + Facilitation sparse 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5479
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4746

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.93
## Success rate for facilitation: 0
data<-sim_data$FacSparseSp10

#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam10fs.rda",interact=fac_inter[[8]])
gjam_mod<-load_gjam(data,name="./gjam_models/gjam10fs.rda", interact=fac_inter[[8]])

#gjfd5<-load_object("./gjam_models/gjam10fs.rda")

g_metric_facSparse10<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[8]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse10$success_env
# g_metric_facSparse10$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 0.326
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacSparseSp10

hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10fs.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[8]])
h_metric_facSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse10$success_env
# h_metric_facSparse10$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse10$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[8]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[8]])

Environmental filtering + Facilitation sparse 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5653
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4826

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.978
## Success rate for facilitation: 0

Gjam

data<-sim_data$FacSparseSp20

#fit_gjam(data,10000,1500,"./gjam_models/gjam20fs.rda",interact=fac_inter[[9]])
gjam_mod <-load_gjam(data,10000,1500,name="./gjam_models/gjam20fs.rda", interact=fac_inter[[9]])

#gjfd5<-load_object("./gjam_models/gjam20fd.rda")

g_metric_facSparse20<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[9]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse20$success_env
# g_metric_facSparse20$success_fac

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 0.446
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0.25
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

HMSC

data<-sim_data$FacSparseSp20
hm_mod<-fit_hmsc(data,nsamples=3000, nchains=2,"Load",name="./HMmodels/hm20fs.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = fac_inter[[9]])
h_metric_facSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[9]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse20$success_env
# h_metric_facSparse20$success_fac
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse20$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[9]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[9]])

Environmental filtering + Competition dense 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## Successful convergence based on Rhat values (all < 1.1).
## Maximum Rhat value for Beta: 1.0324
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.0313

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 1
## [1] 1
## Success rate for non-interacting: 1
## Success rate for facilitation: 0

Gjam

data<-sim_data$CompDenseSp5

#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam5cmpd.rda",interact=comp_inter[[10]])

gjam_mod<-load_gjam(data,10000,1000,name="./gjam_models/gjam5cmpd.rda", interact=(-1)*comp_inter[[10]])

#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")

g_metric_compDense5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[10]], only_env = F)
# cat("Success rate ")
# g_metric_compDense5$success_env
# g_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 0.625
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 
data<-sim_data$CompDenseSp5

hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hm5cmpd.rda" )
#hm_mod<-fit_hmsc(data,"Fit",nsamples=3000, nchains=2 )

hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[10]])

h_metric_compDense5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[10]],only_env = F)
# cat("Success rate ")
# h_metric_compDense5$success_env
# h_metric_compDense5$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_compDense5$success_comp,3)))
## Success rate for facilitation: 1
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[10]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[10]])

Environmental filtering + Competition dense 10 species

## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1845
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1771

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.875
## Success rate for facilitation: 1

Gjam

data<-sim_data$CompDenseSp10

#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10cmpd.rda",interact=comp_inter[[11]])

gjam_mod<-load_gjam(data,5000,500,name="./gjam_models/gjam10cmpd.rda", interact= (-1)*comp_inter[[11]])

#gjfd5<-load_object("./gjam_models/gjam10cmpd.rda")

g_metric_compDense10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[11]], only_env = F)
# cat("Success rate ")
# g_metric_compDense10$success_env
# g_metric_compDense10$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.325
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$CompDenseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm10cmpd.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod,nsamples=5000, nchains=2, interact =  (-1)*comp_inter[[11]])

h_metric_compDense10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[11]],only_env = F)
# cat("Success rate ")
# h_metric_compDense10$success_env
# h_metric_compDense10$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense10$success_env,3)))
## Success rate for non-interacting: 0.875
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense10$success_comp,3)))
## Success rate for competition: 1
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[11]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[11]])

Environmental filtering + Competition dense 20 species

## Summary for model '/tmp/RtmpKix4lZ/file40e96843124a'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2061.424 minutes at time 2019-04-16 03:39:17.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.9065
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.6944

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.956
## Success rate for competition: 0.7

Gjam

data<-sim_data$CompDenseSp20

#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam20cmpd.rda",interact= (-1)*comp_inter[[12]])

gjam_mod<-load_gjam(data,10000,1000,name="./gjam_models/gjam20cmpd.rda", interact= (-1)*comp_inter[[12]])

#gjfd5<-load_object("./gjam_models/gjam20cmpd.rda")

g_metric_compDense20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[12]], only_env = F)
# cat("Success rate ")
# g_metric_compDense20$success_env
# g_metric_compDense20$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.394
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.9
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$CompDenseSp20

hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm20cmpd.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=5000, nchains=2,interact =  (-1)*comp_inter[[12]])

h_metric_compDense20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[12]],only_env = F)
# cat("Success rate ")
# h_metric_compDense20$success_env
# h_metric_compDense20$success_comp
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense20$success_env,3)))
## Success rate for non-interacting: 0.967
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense20$success_comp,3)))
## Success rate for competition: 0.4
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[12]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[12]])

Environmental filtering + Competition sparse 5 species

## Summary for model '/tmp/RtmpKix4lZ/file40e97e225117'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.251 minutes at time 2019-04-17 14:00:46.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1876
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2038

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 1
## Success rate for competition: 1

Gjam

data<-sim_data$CompSparseSp5


#Rho_sign<-fit_gjam(data,5000,2000,"./gjam_models/gjam5cmps.rda",interact= (-1)*comp_inter[[13]])
gjam_mod<-load_gjam(data,5000,2000,name="./gjam_models/gjam5cmps.rda", interact= (-1)*comp_inter[[13]])

#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")

g_metric_compSparse5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[13]], only_env = F)
cat("Success rate ")
## Success rate
g_metric_compSparse5$success_env
## [1] 0.5555556
g_metric_compSparse5$success_comp
## [1] 1
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse5$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$CompSparseSp5

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm5cmps.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[13]])

h_metric_compSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[13]],only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse5$success_comp,3)))
## Success rate for competition: 1
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[13]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[13]])

Environmental filtering + Competition sparse 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e978d641ed'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 685.719 minutes at time 2019-04-17 14:30:01.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.116
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1753

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate
## [1] 0.9534884
## [1] 1
## Success rate for non-interacting: 0.953
## Success rate for competition: 1

Gjam

data<-sim_data$CompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam10cmps.rda",interact= (-1)*comp_inter[[14]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjam10cmps.rda", interact= (-1)*comp_inter[[14]])

g_metric_compSparse10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[14]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse10$success_env
# g_metric_compSparse10$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse10$success_env,3)))
## Success rate for non-interacting: 0.326
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse10$success_comp,3)))
## Success rate for competition: 1
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 
data<-sim_data$CompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm10cmps.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[14]])

h_metric_compSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[14]],only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse10$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse10$success_comp,3)))
#mod_list<-list.append(mod_list,hmsc=Rho_sign)

mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[14]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[14]])

Environmental filtering + Competition sparse 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9379ee963'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2084.707 minutes at time 2019-04-18 01:55:46.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2652
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3883

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.968
## Success rate for competition: 0

Gjam

data<-sim_data$CompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam20cmps.rda",interact= (-1)*comp_inter[[15]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjam20cmps.rda", interact= (-1)*comp_inter[[15]])

g_metric_compSparse20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[15]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse20$success_env
# g_metric_compSparse20$success_comp

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 0.269
cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0.75
#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$CompSparseSp20

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm20cmps.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[15]])

h_metric_compSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[15]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse20$success_comp,3)))
## Success rate for competition: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)

mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[15]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[15]])

Environmental filtering + Competition +Facilitation 5 dense species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e93b49ad29'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 28.267 minutes at time 2019-04-19 12:40:31.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1667
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2086

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.667
## Success rate for competition: 1
## Success rate for facilitation: 0.5

Gjam

data<-sim_data$FacCompDenseSp5
#Rho_sign<-fit_gjam(data,30000,15000,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
gjam_mod<-load_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])

#Rho_sign3<-fit_gjam(data,5000,1000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign4<-fit_gjam(data,5000,2500,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign5<-fit_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmp_facd.rda", (-1)*comp_inter[[16]]+fac_inter[[16]])

g_metric_FacCompDense5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[16]], fac=fac_inter[[16]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense5$success_fac,3)))
## Success rate for facilitation: 1
# par(mfrow=c(2,3))
#corrplot((Rho_sign+Rho_sign2)/2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign3, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign4, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign5, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")


#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacCompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hmcmp_facd5.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact =  (-1)*comp_inter[[16]]+ fac_inter[[16]])

h_metric_FacCompDense5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[16]],fac= fac_inter[[16]],only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense5$success_env,3)))
## Success rate for non-interacting: 0.833
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense5$success_comp,3)))
## Success rate for competition: 0.5
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense5$success_fac,3)))
## Success rate for facilitation: 0.5
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])

Environmental filtering + Competition +Facilitation 10 dense species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e93b173410'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 688.628 minutes at time 2019-04-19 13:08:47.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.3458
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3195

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.914
## Success rate for competition: 0
## Success rate for facilitation: 0

Gjam

data<-sim_data$FacCompDenseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd10.rda",interact=  (-1)*comp_inter[[17]]+fac_inter[[17]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facd10.rda", interact= (-1)*comp_inter[[17]]+fac_inter[[17]])

g_metric_FacCompDense10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[17]], fac=fac_inter[[17]], only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense10$success_env,3)))
## Success rate for non-interacting: 0.286
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense10$success_comp,3)))
## Success rate for competition: 0.6
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense10$success_fac,3)))
## Success rate for facilitation: 1
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd10.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[17]] + fac_inter[[17]])

h_metric_FacCompDense10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[17]],fac=fac_inter[[17]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense10$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])

Environmental filtering + Competition +Facilitation 20 dense species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e95eadc1e3'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2073.109 minutes at time 2019-04-20 00:37:26.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.3008
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2755

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.965
## Success rate for competition: 0.333
## Success rate for facilitation: 0.333

Gjam

data<-sim_data$FacCompDenseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd20.rda",interact=  (-1)*comp_inter[[18]]+fac_inter[[18]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facd20.rda", interact= (-1)*comp_inter[[18]]+fac_inter[[18]])

g_metric_FacCompDense20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[18]], fac=fac_inter[[18]], only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense20$success_env,3)))
## Success rate for non-interacting: 0.366
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense20$success_comp,3)))
## Success rate for competition: 0.667
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense20$success_fac,3)))
## Success rate for facilitation: 0.429
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacCompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd20.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[18]] + fac_inter[[18]])
h_metric_FacCompDense20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[18]],fac=fac_inter[[18]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense20$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense20$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense20$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])

Environmental filtering + Competition +Facilitation 5 sparse species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9896a5fb'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.285 minutes at time 2019-04-21 11:10:35.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.3771
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3486

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.875
## Success rate for competition: 1
## Success rate for facilitation: 1

Gjam

data<-sim_data$FacCompSparseSp5
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs5.rda",interact= (-1)*comp_inter[[19]]+ fac_inter[[19]] )
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs5.rda", interact=  (-1)*comp_inter[[19]]+ fac_inter[[19]])

g_metric_FacCompSparse5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[15]], only_env = F)

cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse5$success_env,3)))
## Success rate for non-interacting: 0.556
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse5$success_comp,3)))
## Success rate for competition: 1
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacCompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs5.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[19]]+ fac_inter[[19]])
 
h_metric_FacCompSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[19]],fac= fac_inter[[19]] ,only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse5$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse5$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse5$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])

Environmental filtering + Competition +Facilitation 10 sparse species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e95f910eba'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 679.081 minutes at time 2019-04-21 11:39:53.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 2.1146
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 3.3048

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.951
## Success rate for competition: 0
## Success rate for facilitation: 0

Gjam

data<-sim_data$FacCompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs10.rda",interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs10.rda", interact=  (-1)*comp_inter[[20]] +fac_inter[[20]])

g_metric_FacCompSparse10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]], only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse10$success_env,3)))
## Success rate for non-interacting: 0.341
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse10$success_comp,3)))
## Success rate for competition: 1
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparse10$success_fac,3)))
## Success rate for facilitation: 0
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacCompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs10.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[20]] +fac_inter[[20]])

h_metric_FacCompSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse10$success_env,3)))
## Success rate for non-interacting: 1
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse10$success_comp,3)))
## Success rate for competition: 0
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse10$success_fac,3)))
## Success rate for facilitation: 0
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])

Environmental filtering + Competition +Facilitation 20 sparse species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e978c52358'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2083.296 minutes at time 2019-04-21 22:58:59.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2707
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3875

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
## Success rate for non-interacting: 0.984
## Success rate for competition: 0
## Success rate for facilitation: 0

Gjam

data<-sim_data$FacCompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs20.rda",interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs20.rda", interact= (-1)*comp_inter[[21]]+fac_inter[[21]])

g_metric_FacCompSparse20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[21]], fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse20$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse20$success_comp,3)))
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparse20$success_fac,3)))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign) 
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) 

HMSC

data<-sim_data$FacCompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs20.rda" )
hm_conv(hm_mod)

hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[21]] +fac_inter[[21]])
h_metric_FacCompSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[21]],fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse20$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse20$success_comp,3)))
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse20$success_fac,3)))
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)

ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])

ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])

Graph from paper

Visualization for true interactions

Facilitation

Competition Dense

Competition Dense

Competition + Facilitation Dense

Competition + Facilitation Sparse